{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 基于梯度算法的脉冲优化\n",
    "\n",
    "*版权所有 (c) 2021 百度量子计算研究所，保留所有权利。*"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 内容概要\n",
    "本教程介绍如何通过梯度上升脉冲工程（GRadient Ascent Pulse Engineering，GRAPE）算法产生高保真度的单量子比特门脉冲。本教程的大纲如下:\n",
    "\n",
    "- 背景介绍\n",
    "- 准备工作\n",
    "- 构建哈密顿量\n",
    "- 生成目标门的 GRAPE 优化脉冲\n",
    "- 总结"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 背景介绍\n",
    "\n",
    "**GRAPE 算法的原理**\n",
    "\n",
    "梯度算法的主要目的是通过迭代找到目标函数的最值，或者收敛目标函数到最值。若我们想寻找最大值，则称之为梯度上升算法：\n",
    "$$\n",
    "\\Theta_n = \\Theta_{n-1} + k\\bigtriangledown{J(\\Theta)}, \\tag{1}\n",
    "$$\n",
    "\n",
    "梯度上升算法和“爬山”非常相似，若我们想要爬到“山顶”，需要知道每次迈步的方向和步长 $k\\bigtriangledown{J(\\Theta)}$。根据上式，每次“迈步”后的位置 $\\Theta_n$ 由原位置 $\\Theta_{n-1}$、目标函数 $J(\\Theta)$ 在原位置的导数和系数 $k$ 的乘积决定。\n",
    "\n",
    "GRAPE 算法 \\[1\\] 是量子优化控制领域中的经典方法，最初在核磁共振平台中得到应用，后来扩展到了其它平台中。其核心思想是，通过对量子门的执行时间进行切片，并假设每片内的脉冲为常量，随后将每一切片中脉冲的强度作为优化参数，并将薛定谔方程作为优化的限制条件进行梯度上升优化求解。具体来说，量子态的动力学演化规律遵循薛定谔方程，因此我们可以通过求解海森堡绘景 (Heisenberg picture) 中的薛定谔方程来得到时间区间 $t \\in \\left[0, T\\right]$ 内系统演化所对应的演化算符 $U (t)$：\n",
    "$$\n",
    "i\\hbar \\frac{\\partial U (t)}{\\partial t} = \\hat{H}(t)U (t),\n",
    "\\tag{2}\n",
    "$$\n",
    "\n",
    "其中 $\\hbar$ 是约化普朗克常数。上式是一个微分方程，在计算机上通常使用数值求解的方式近似求解。当哈密顿量不含时，即 $\\hat{H}(t)$ 不随时间 $t$ 变化时，可以使用矩阵指数求解系统的演化算符（从这里开始，我们取 $\\hbar=1$）：\n",
    "\n",
    "$$\n",
    "U(T) = \\exp(-i\\hat{H}T).\n",
    "\\tag{3}\n",
    "$$\n",
    "\n",
    "但对于含时（并且各个时刻不对易）的哈密顿量而言，上述方法不再适用。一种有效的方法是对整个演化过程在时间上切片，并假设在每一个切片 $\\left[t_j, t_{j+1}\\right]$ 内， 哈密顿量 $H_{j}$ 是不含时的，其中下标 $j$ 为切片的顺序编号。因此每一个切片对应的演化算符 $U_{j}$ 可以使用式 $(2)$ 进行求解。最终整个演化过程可以写为：\n",
    "\n",
    "$$\n",
    "U = U_{N}U_{N-1}\\cdots{U}_{1},\n",
    "\\tag{4}\n",
    "$$\n",
    "\n",
    "其中 $N$ 为切片的数量。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**GRAPE 算法的实现**\n",
    "\n",
    "下面我们介绍以未归一的门保真度的平方作为目标函数 $J$ 的 GRAPE 实现过程：\n",
    "\n",
    "$$\n",
    "J = \\left|{\\rm Tr}\\left[(U^{\\dagger}_{\\rm target}U(T)\\right]\\right|^2 = \\left|{\\rm Tr}\\left[P^{\\dagger}_{j}X_{j}\\right]\\right|^2,  \\ \\forall\\ 0<j<N, \\tag{5}\n",
    "$$\n",
    "\n",
    "其中 $U_{\\rm target}$ 是目标单量子比特门，$X_{j} = U_{j}\\cdots{U}_{1}$ 是中间传播算子，$P_{j} = U^{\\dagger}_{j+1}\\cdots{U}^{\\dagger}_{N}U_{\\rm target}$ 是中间反向传播算子。\n",
    "\n",
    "连续的脉冲波形 $u_i(t)$ 可以切片后由离散的 $u_i(j)$ 表示，其中 $i$ 指标表示不同控制项（例如当 $i=x$ 时 $\\hat{H}_{i}=\\hat{H}_{x}$）上的脉冲。对目标函数 $J$ 求偏导可得 \\[2\\]：\n",
    "\n",
    "$$\n",
    "\\dfrac{\\partial J}{\\partial u_i(j)} = {\\rm Re}\\{-2 i \\Delta t \\left[{\\rm Tr}(P^{\\dagger}_{j}\\hat{H}_iX_{j}){\\rm Tr}(P^{\\dagger}_{j}X_{j})\\right]\\},\n",
    "\\tag{6}\n",
    "$$\n",
    "\n",
    "根据基于梯度的算法的原理，当设置了学习率 $k$ 后，每次迭代优化时算法可对脉冲进行“整形”，最终将优化出一个不规则的脉冲：\n",
    "$$\n",
    "u_i(t) \\mapsto u_i(t) + k\\dfrac{\\partial J}{\\partial u_i(j)},\\tag{7}\n",
    "$$\n",
    "\n",
    "![GRAPE](figures/GRAPE1.png)\n",
    "\n",
    "如上图所示，GRAPE 的每一次迭代相当于对脉冲波形进行了一次“整形”。\n",
    "\n",
    "下面我们将介绍如何利用量脉实现 GRAPE 这一过程。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 准备工作\n",
    "\n",
    "您可以按照本教程运行下面的程序。在成功安装量脉后，需要从量脉（Quanlse）和其它常用 Python 库导入以下包："
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Import required packages\n",
    "from Quanlse.remoteOptimizer import remoteOptimize1QubitGRAPE\n",
    "from Quanlse.QOperator import number, driveX, driveY, duff\n",
    "from Quanlse.QHamiltonian import QHamiltonian as QHam\n",
    "from Quanlse.QOperation import FixedGate\n",
    "from Quanlse.Utils.Functions import project\n",
    "from Quanlse.QWaveform import gaussian\n",
    "\n",
    "import numpy as np\n",
    "from numpy import dot, round"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "首先，我们设置哈密顿量的基本参数，在这里我们考虑一个二能级的单量子比特系统："
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Sampling period\n",
    "dt = 1.0\n",
    "\n",
    "# System energy level\n",
    "level = 2\n",
    "\n",
    "# Duration of the gate (ns)\n",
    "tg = 40"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 构造哈密顿量\n",
    "\n",
    "接下来，我们定义系统的哈密顿量。在旋转坐标系（Rotating Frame）转换后，描述该量子系统的哈密顿量可以写为：\n",
    "\n",
    "$$\n",
    "\\hat{H} = \\frac{1}{2} \\Omega^x(t) (\\hat{a}+\\hat{a}^{\\dagger}) + i \\frac{1}{2} \\Omega^y(t) (\\hat{a}-\\hat{a}^{\\dagger}) ,\n",
    "$$\n",
    "\n",
    "其中 $\\Omega^x(t)$ 是 X 通道的微波脉冲的强度；$\\Omega^y(t)$ 是 Y 通道的微波脉冲强度。在量脉中对哈密顿量的构建如下所示："
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Initialize the Hamiltonian\n",
    "ham = QHam(subSysNum=1, sysLevel=level, dt=dt)\n",
    "\n",
    "# Add the anharmonicity term\n",
    "alphaQ = - 0.22 * (2 * np.pi)\n",
    "ham.addDrift(duff, 0, coef=alphaQ)\n",
    "\n",
    "# Add the control terms\n",
    "ham.appendWave(driveX, 0, waves=gaussian(tg, a=0.3, tau=tg / 2, sigma=tg / 8))\n",
    "ham.appendWave(driveY, 0, waves=gaussian(tg, a=0.3, tau=tg / 2, sigma=tg / 8))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "使用量脉云服务之前，我们需要登录 http://quantum-hub.baidu.com 获取 token 来访问云端。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [],
   "source": [
    "# Import Define class and set the token\n",
    "# Please visit http://quantum-hub.baidu.com\n",
    "from Quanlse import Define\n",
    "Define.hubToken = ''"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "pycharm": {
     "name": "#%% md\n"
    }
   },
   "source": [
    "## 生成目标门的 GRAPE 优化脉冲\n",
    "\n",
    "接下来我们调用量脉优化器中的 `remoteOptimize1QubitGRAPE()` 来得到 GRAPE 算法优化出的脉冲与门保真度。在 `remoteOptimize1QubitGRAPE()` 中，用户可以自定义门的持续时间 `tg`（默认为 20 ns）、最大迭代次数 `iterate`（默认为 150）以及优化通道 `xyzPulses`（默认为 X, Y 通道，即 \\[1, 1, 0\\]）。切片数取决于 $tg / dt$。在这里，我们希望实现 H 门：\n",
    "$$\n",
    "U_{\\rm target} = \n",
    "\\dfrac{1}{\\sqrt{2}}\\begin{bmatrix} \n",
    "1 & 1  \\\\ 1 & -1 \n",
    "\\end{bmatrix} .\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "pycharm": {
     "name": "#%%\n"
    },
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "# Define the target unitary evolution\n",
    "uGoal = FixedGate.H.getMatrix()\n",
    "\n",
    "# Run the optimization\n",
    "job, infid = remoteOptimize1QubitGRAPE(ham, uGoal, tg=tg, iterate=50, xyzPulses=None)\n",
    "\n",
    "# Print infidelity and the waveforms\n",
    "print(f\"minimum infidelity with GRAPE: {infid}\")\n",
    "ham.plot(color = ['blue', 'mint'], dark=True)\n",
    "\n",
    "# Print the evolution\n",
    "result = ham.simulate()\n",
    "print(\"The evolution U:\\n\", round(result.result[0][\"unitary\"], 2))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "经过 GRAPE 算法优化出的量子门与目标量子门通常相差一个全局相位 $e^{i\\phi}$，从作用效果出发，这个全局相位是可以忽略的。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 总结\n",
    "\n",
    "本教程介绍了使用 GRAPE 算法生成单量子比特门脉冲。我们可以看到，GRAPE 算法优化出来的脉冲波形没有规则的形状。作为最常见的优化算法，GRAPE 算法在合理的初始值下收敛速度很快。GRAPE算法的效果与切片数量有关。我们鼓励用户尝试不同于本教程的参数值与优化器 `remoteOptimizer1Qubit` 多作对比以获得最佳结果。\n",
    "\n",
    "用户可以点击这个链接 [tutorial-GRAPE.ipynb](https://github.com/baidu/Quanlse/blob/main/Tutorial/CN/tutorial-GRAPE-cn.ipynb) 跳转到此 Jupyter Notebook 文档相应的 GitHub 页面来获取相关代码并进一步探索 GRAPE 算法。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 参考资料\n",
    "\n",
    "\\[1\\] [Wilhelm, Frank K., et al. \"An introduction into optimal control for quantum technologies.\" *arXiv preprint arXiv*:2003.10132 (2020).](https://arxiv.org/abs/2003.10132v1)\n",
    "\n",
    "\\[2\\] [Khaneja, Navin, et al. \"Optimal control of coupled spin dynamics: design of NMR pulse sequences by gradient ascent algorithms.\" *Journal of magnetic resonance* 172.2 (2005): 296-305.](https://doi.org/10.1016/j.jmr.2004.11.004)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
